## loading libraries
library(ggplot2)
library(ggpubr)
## Loading required package: magrittr
library(Seurat)
library(cowplot)
##
## ********************************************************
## Note: As of version 1.0.0, cowplot does not change the
## default ggplot2 theme anymore. To recover the previous
## behavior, execute:
## theme_set(theme_cowplot())
## ********************************************************
##
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggpubr':
##
## get_legend
library(readxl)
library(corto)
library(Matrix)
library(EnhancedVolcano)
## Loading required package: ggrepel
## Initilization
set.seed(0)
label_size=8
point_size=1
folder="~/corona_project/plots_intigrate/"
cell_type_genes=list("Bowman’s glands"=c("SOX9", "SOX10", "GPX3"),
"olfactory HBCs"=c("TP63", "KRT5", "CXCL14", "SOX2", "MEG3"),
"olfactory ensheathing glia"=c("S100B", "PLP1", "PMP2","MPZ", "ALX3"),
"olfactory microvillar cells"=c("ASCL3", "CFTR", "HEPACAM2", "FOXL1"),
"immature neurons"=c("GNG8", "OLIG2", "EBF2", "LHX2", "CBX8"),
"mature neurons"=c("GNG13", "EBF2", "CBX8", "RTP1"),
"GBCs"=c("HES6", "ASCL1", "CXCR4", "SOX2", "EZH2", "NEUROD1", "NEUROG1"),
"sustentacular cells"=c("CYP2A13", "CYP2J2", "GPX6", "ERMN", "SOX2"))
cell_type_genes_temp=list("Bowman’s glands"=c("SOX9", "GPX3"),
"olfactory HBCs"=c("CXCL14", "MEG3"),
"olfactory ensheathing glia"=c("S100B", "PLP1"),
"olfactory microvillar cells"=c( "ASCL3", "HEPACAM2"),
"immature neurons"=c("GNG8","OLIG2"),
"mature neurons"=c("GNG13","RTP1"),
"GBCs"=c("CXCR4", "SOX2"),
"sustentacular cells"=c("CYP2A13", "ERMN"))
marker_genes=c("SOX9", "SOX10", "GPX3",
"TP63", "KRT5", "CXCL14", "SOX2", "MEG3",
"S100B", "PLP1", "PMP2","MPZ", "ALX3",
"ASCL3", "CFTR", "HEPACAM2", "FOXL1",
"GNG8", "OLIG2", "EBF2", "LHX2", "CBX8",
"GNG13", "EBF2", "CBX8", "RTP1",
"HES6", "ASCL1", "CXCR4", "SOX2", "EZH2", "NEUROD1", "NEUROG1",
"CYP2A13", "CYP2J2", "GPX6", "ERMN", "SOX2")
marker_genes_temp=c("SOX9", "GPX3","CXCL14", "MEG3","S100B", "PLP1", "ASCL3", "HEPACAM2",
"GNG8","OLIG2","GNG13","RTP1","CXCR4", "SOX2","CYP2A13", "ERMN")
cell_type_names=c("Bowman’s glands",
"olfactory HBCs",
"olfactory ensheathing glia",
"olfactory microvillar cells",
"immature neurons",
"mature neurons",
"GBCs",
"sustentacular cells")
## Laoding paitent P2 and P3 matrices and formation of seurat object
batch_list=list("P2","P3")
batch_data_list=list("P2"=1,"P3"=1)
for( i in 1:length(batch_list))
{
print(batch_list[[i]])
s_object=Read10X(paste("~/corona_project/Input_files/",batch_list[[i]],sep=""))
s_object=CreateSeuratObject(counts =s_object, min.cells = 0, min.features = 400, project = "P23")
s_object[["percent.mt"]] <- PercentageFeatureSet(s_object, pattern = "^MT-")
s_object <- subset(s_object, subset = nFeature_RNA >100 & nFeature_RNA <8000 & percent.mt <10)
s_object@meta.data[, "run"] <- batch_list[i]
s_object=NormalizeData(s_object)
batch_data_list[[i]]=FindVariableFeatures(s_object, selection.method = "vst", nfeatures =5000)
}
## [1] "P2"
## [1] "P3"
batch_data_list
## $P2
## An object of class Seurat
## 33538 features across 10881 samples within 1 assay
## Active assay: RNA (33538 features, 5000 variable features)
##
## $P3
## An object of class Seurat
## 33538 features across 5415 samples within 1 assay
## Active assay: RNA (33538 features, 5000 variable features)
saveRDS(batch_data_list,"~/corona_project/batch_data_list_5000_23.rds")
## Intigration and batch effect removing using Seurat
set.seed(1)
run.anchors <- FindIntegrationAnchors(object.list = batch_data_list, dims = 1:30,anchor.features = 5000)
## Computing 5000 integration features
## Scaling features for provided objects
## Finding all pairwise anchors
## Running CCA
## Merging objects
## Finding neighborhoods
## Finding anchors
## Found 11489 anchors
## Filtering anchors
## Retained 8787 anchors
## Extracting within-dataset neighbors
a=run.anchors
run.combined <- IntegrateData(anchorset = run.anchors, dims = 1:30)
## Merging dataset 2 into 1
## Extracting anchors for merged samples
## Finding integration vectors
## Finding integration vector weights
## Integrating data
## Warning: Adding a command log without an assay associated with it
b=run.combined
sce_3_1_1_before=run.combined
saveRDS(sce_3_1_1_before,"~/corona_project/sce_3_1_1_before_5000_23.rds")
## Dim reduction and clustering
DefaultAssay(run.combined) <- "integrated"
run.combined <- ScaleData(run.combined, verbose = FALSE)
run.combined <- RunPCA(run.combined, npcs = 30, verbose = FALSE)
run.combined <- RunUMAP(run.combined, reduction = "pca", dims = 1:30)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 23:48:30 UMAP embedding parameters a = 0.9922 b = 1.112
## 23:48:30 Read 16296 rows and found 30 numeric columns
## 23:48:30 Using Annoy for neighbor search, n_neighbors = 30
## 23:48:30 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 23:48:33 Writing NN index file to temp file /tmp/RtmpQ1SbcN/file18d38f17eda
## 23:48:33 Searching Annoy index using 1 thread, search_k = 3000
## 23:48:37 Annoy recall = 100%
## 23:48:37 Commencing smooth kNN distance calibration using 1 thread
## 23:48:38 Initializing from normalized Laplacian + noise
## 23:48:39 Commencing optimization for 200 epochs, with 707164 positive edges
## 23:48:44 Optimization finished
run.combined <- FindNeighbors(run.combined, reduction = "pca", dims = 1:30)
## Computing nearest neighbor graph
## Computing SNN
run.combined <- FindClusters(run.combined, resolution = 0.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 16296
## Number of edges: 654023
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9526
## Number of communities: 22
## Elapsed time: 2 seconds
## Cell-types annotation
for(i in 0:as.integer(length(marker_genes)/4))
{
a=(i*4+1)
b=((i+1)*4)
if(i==as.integer(length(marker_genes)/4))
b=length(marker_genes)
if(a>length(marker_genes))
next
print(i)
print(VlnPlot(run.combined, features = marker_genes[a:b],ncol=2))
print(FeaturePlot(run.combined, features = marker_genes[a:b],cols=c("lightgrey", "red"),label = TRUE,pt.size=2,order = TRUE))
}
## [1] 0


## [1] 1


## [1] 2


## [1] 3

## [1] 4
## Warning: Could not find FOXL1 in the default search locations, found in RNA
## assay instead

## Warning: Could not find FOXL1 in the default search locations, found in RNA
## assay instead

## [1] 5
## Warning: Could not find CBX8 in the default search locations, found in RNA assay
## instead

## Warning: Could not find CBX8 in the default search locations, found in RNA assay
## instead

## [1] 6
## Warning: Could not find CBX8 in the default search locations, found in RNA assay
## instead

## Warning: Could not find CBX8 in the default search locations, found in RNA assay
## instead


## [1] 7


## [1] 8


## [1] 9


run.combined<-RenameIdents(run.combined,`13`="Immature neurons", `14`="Mature neurons", `3`="Olfactory HBCs",`19`="Olfactory microvillar cells",
`5`="Bowman's gland", `16`="Olfactory ensheathing glia", `9`="Sustentacular cells",`20`="GBCs" )
run.combined=subset(run.combined, idents = c("Immature neurons","Mature neurons","Olfactory HBCs","Olfactory microvillar cells","Bowman's gland",
"Olfactory ensheathing glia","Sustentacular cells","GBCs"))
sce_3_1_1_after=run.combined
saveRDS(sce_3_1_1_after,"~/corona_project/sce_3_1_1_after_5000_23.rds")
## Ploting heat map for selected gene markers on all clusters
sce_3_1_1_after_5000_23 <- readRDS("~/corona_project/sce_3_1_1_after_5000_23.rds")
cluster.averages <- AverageExpression(sce_3_1_1_after_5000_23, return.seurat = TRUE)
## Finished averaging RNA for cluster Immature neurons
## Finished averaging RNA for cluster Mature neurons
## Finished averaging RNA for cluster Olfactory HBCs
## Finished averaging RNA for cluster Olfactory microvillar cells
## Finished averaging RNA for cluster Bowman's gland
## Finished averaging RNA for cluster Olfactory ensheathing glia
## Finished averaging RNA for cluster Sustentacular cells
## Finished averaging RNA for cluster GBCs
## Finished averaging integrated for cluster Immature neurons
## Finished averaging integrated for cluster Mature neurons
## Finished averaging integrated for cluster Olfactory HBCs
## Finished averaging integrated for cluster Olfactory microvillar cells
## Finished averaging integrated for cluster Bowman's gland
## Finished averaging integrated for cluster Olfactory ensheathing glia
## Finished averaging integrated for cluster Sustentacular cells
## Finished averaging integrated for cluster GBCs
## Centering and scaling data matrix
DoHeatmap(sce_3_1_1_after_5000_23,features =c(marker_genes_temp),size = 3, draw.lines = FALSE )

DoHeatmap(cluster.averages,features =c(marker_genes_temp),size = 3, draw.lines = FALSE )

## Vln plots with all selected gene markers
plot=list()
for(i in 1:length(marker_genes_temp))
{
print(i)
if(i %in% c(1,5,9))
plot[[i]]<-print(VlnPlot(sce_3_1_1_after_5000_23, features = marker_genes_temp[[i]])) +
NoLegend() + theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.line.x = element_blank())
if(i %in% c(14,15,16))
plot[[i]]<-print(VlnPlot(sce_3_1_1_after_5000_23, features = marker_genes_temp[[i]])) +
NoLegend() + theme(axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
axis.line.y = element_blank())
if(i %in% c(13))
plot[[i]]<-print(VlnPlot(sce_3_1_1_after_5000_23, features = marker_genes_temp[[i]])) +
NoLegend()
if(i %in% c(2,3,4,6,7,8,10,11,12))
plot[[i]]<-print(VlnPlot(sce_3_1_1_after_5000_23, features = marker_genes_temp[[i]])) +
NoLegend()+ theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.line.x = element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
axis.line.y = element_blank())
}
## [1] 1

## [1] 2

## [1] 3

## [1] 4

## [1] 5

## [1] 6

## [1] 7

## [1] 8

## [1] 9

## [1] 10

## [1] 11

## [1] 12

## [1] 13

## [1] 14

## [1] 15

## [1] 16

cowplot::plot_grid(plotlist = plot,nrow=4,rel_heights = c(1,1,1,2))

## Feature plots with all selected gene markers
plot=list()
plot[[1]]=DimPlot(sce_3_1_1_after_5000_23, reduction = "umap",pt.size = point_size,label=TRUE)+
labs(title = "") + NoAxes() + NoLegend()
for(i in 1:length(marker_genes_temp))
{
print(i)
plot[[i+1]]<-print(FeaturePlot(sce_3_1_1_after_5000_23, features = marker_genes_temp[[i]],cols=c("lightgrey", "red"),
pt.size = point_size,label.size = label_size,order = TRUE)) +
NoLegend() + NoAxes()
}
## [1] 1

## [1] 2

## [1] 3

## [1] 4

## [1] 5

## [1] 6

## [1] 7

## [1] 8

## [1] 9

## [1] 10

## [1] 11

## [1] 12

## [1] 13

## [1] 14

## [1] 15

## [1] 16

AT=cowplot::plot_grid(plotlist = plot[c(3,4,5)],ncol=3)
AB=cowplot::plot_grid(plotlist = plot[c(12,13,11)],ncol=3)
AR=cowplot::plot_grid(plotlist = plot[c(6,7,8,9,10)],nrow=5)
AL=cowplot::plot_grid(plotlist = plot[c(2,14,15,16,17)],nrow=5)
AC=cowplot::plot_grid(plotlist = plot[c(1)],nrow=1)
ATCB=cowplot::plot_grid(plotlist = list(AT,AC,AB),nrow=3,rel_heights = c(1,4,1))
cowplot::plot_grid(plotlist = list(AL,ATCB,AR),ncol=3,rel_widths = c(1,2,1))

### Dim plot to show all cell types with different colors
DimPlot(sce_3_1_1_after_5000_23, reduction = "umap",pt.size = point_size,label=TRUE)+ labs(title = "") + NoLegend()

DimPlot(sce_3_1_1_after_5000_23, reduction = "umap",pt.size = point_size,label=FALSE)+ labs(title = "")

## FeaturePlot to show expression of gene markers in respective cell markers
### FeaturePlot to show the expression of ACE2 in each cell type
FeaturePlot(run.combined, features = "ACE2",combine = FALSE,cols=c("lightgrey", "red"),pt.size = point_size,label.size = label_size,order =TRUE)
## Warning: Could not find ACE2 in the default search locations, found in RNA assay
## instead
## [[1]]

### FeaturePlot to show the expression of TMPRSS2 in each cell type
FeaturePlot(run.combined, features = "TMPRSS2",combine = FALSE,cols=c("lightgrey", "green"),pt.size = point_size,label.size = label_size,order =TRUE)
## [[1]]

### FeaturePlot to show the expression of CTSL in each cell type
FeaturePlot(run.combined, features = "CTSL",combine = FALSE,cols=c("lightgrey", "blue"),pt.size = point_size,label.size = label_size,order =TRUE)
## Warning: Could not find CTSL in the default search locations, found in RNA assay
## instead
## [[1]]

### FeaturePlot to show the expression of BSG in each cell type
FeaturePlot(run.combined, features = "BSG",combine = FALSE,cols=c("lightgrey", "magenta"),pt.size = point_size,label.size = label_size,order =TRUE)
## [[1]]

## Initiallizing a matrix for 24 * 8, 8 for cells type and 24 for different cases as shown below
count_matrix=matrix(0,24,8)
rownames(count_matrix)<-c("ACE2_count","ACE2_percent","TMPRSS2_count","TMPRSS2_percent",
"CTSL_count","CTSL_percent","BSG_count","BSG_percent",
"intersect_of_ACE2_TMPRSS2_count","ACE2_contri_percent","TMPRSS2_contri_percent",
"intersect_of_ACE2_CTSL_count","ACE2_contri_percent","CTSL_contri_percent",
"intersect_of_BSG_TMPRSS2_count","BSG_contri_percent","TMPRSS2_contri_percent",
"intersect_of_BSG_CTSL_count","BSG_contri_percent","CSTL_contri_percent",
"ACE2_TMPRSS2_CTSL","ACE2_TMPRSS2_CTSB",
"BSB_TMPRSS2_CTSL","BSG_TMPRSS2_CTSB")
colnames(count_matrix)<-unique(Idents(run.combined))
### Percentage of ACE2 expressed's cells in each cell type
a=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["ACE2",])>0])])
count_matrix[1,names(a)]=a
print(sum(a))
## [1] 32
a=a*100/sum(a)
count_matrix[2,names(a)]=a
df=data.frame(cell_type=names(a),count_percentage=a)
ggplot(data=df, aes(x=cell_type, y=count_percentage.Freq, fill=cell_type)) +
geom_bar(stat="identity") + theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1),legend.position="none")

dim(run.combined@assays$RNA)
## [1] 33538 3906
### Percentage of TMPRSS2 expressed's cells in each cell type
a=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["TMPRSS2",])>0])])
count_matrix[3,names(a)]=a
print(sum(a))
## [1] 666
a=a*100/sum(a)
count_matrix[4,names(a)]=a
df=data.frame(cell_type=names(a),count_percentage=a)
ggplot(data=df, aes(x=cell_type, y=count_percentage.Freq, fill=cell_type)) +
geom_bar(stat="identity") + theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1),legend.position="none")

dim(run.combined@assays$RNA)
## [1] 33538 3906
### Percentage of CTSL expressed's cells in each cell type
a=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["CTSL",])>0])])
count_matrix[5,names(a)]=a
print(sum(a))
## [1] 1631
a=a*100/sum(a)
count_matrix[6,names(a)]=a
df=data.frame(cell_type=names(a),count_percentage=a)
ggplot(data=df, aes(x=cell_type, y=count_percentage.Freq, fill=cell_type)) +
geom_bar(stat="identity") + theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1),legend.position="none")

dim(run.combined@assays$RNA)
## [1] 33538 3906
### Percentage of BSG expressed's cells in each cell type
a=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["BSG",])>0])])
count_matrix[7,names(a)]=a
print(sum(a))
## [1] 3014
a=a*100/sum(a)
count_matrix[8,names(a)]=a
df=data.frame(cell_type=names(a),count_percentage=a)
ggplot(data=df, aes(x=cell_type, y=count_percentage.Freq, fill=cell_type)) +
geom_bar(stat="identity") + theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1),legend.position="none")

dim(run.combined@assays$RNA)
## [1] 33538 3906
### Bar plot of intersect of both ACE2 and TMPRSS2
col_name_TMPRSS2=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["TMPRSS2",])>0])]
col_name_ACE2=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["ACE2",])>0])]
intersect_names=intersect(names(col_name_ACE2),names(col_name_TMPRSS2))
intersect_cells=table(Seurat::Idents(run.combined)[intersect_names])
count_matrix[9,names(intersect_cells)]=intersect_cells
TMPRSS2_cells=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["TMPRSS2",])>0])])
ACE2_cells=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["ACE2",])>0])])
count_matrix[10,names(ACE2_cells)]=(intersect_cells/ACE2_cells)*100
count_matrix[11,names(TMPRSS2_cells)]=(intersect_cells/TMPRSS2_cells)*100
new_data=data.frame("cells_percenrage"=c(((intersect_cells/TMPRSS2_cells)*100),((intersect_cells/ACE2_cells)*100)),
genes=c(rep("TMPRSS2",length(TMPRSS2_cells)),rep("ACE2",length(ACE2_cells))),
cell_types=c(names(TMPRSS2_cells),names(ACE2_cells)))
ggplot(data=new_data, aes(x=cell_types, y=cells_percenrage, fill=genes)) +
geom_bar(stat="identity") + theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))
## Warning: Removed 4 rows containing missing values (position_stack).

### Bar plot of intersect of both ACE2 and CTSL
col_name_CTSL=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["CTSL",])>0])]
col_name_ACE2=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["ACE2",])>0])]
intersect_names=intersect(names(col_name_ACE2),names(col_name_CTSL))
intersect_cells=table(Seurat::Idents(run.combined)[intersect_names])
count_matrix[12,names(intersect_cells)]=intersect_cells
CTSL_cells=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["CTSL",])>0])])
ACE2_cells=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["ACE2",])>0])])
count_matrix[13,names(ACE2_cells)]=(intersect_cells/ACE2_cells)*100
count_matrix[14,names(CTSL_cells)]=(intersect_cells/CTSL_cells)*100
new_data=data.frame("cells_percenrage"=c(((intersect_cells/CTSL_cells)*100),((intersect_cells/ACE2_cells)*100)),
genes=c(rep("CTSL",length(CTSL_cells)),rep("ACE2",length(ACE2_cells))),
cell_types=c(names(CTSL_cells),names(ACE2_cells)))
ggplot(data=new_data, aes(x=cell_types, y=cells_percenrage, fill=genes)) +
geom_bar(stat="identity") + theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))
## Warning: Removed 4 rows containing missing values (position_stack).

### Bar plot of intersect of both BSG and TMPRSS2
col_name_TMPRSS2=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["TMPRSS2",])>0])]
col_name_BSG=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["BSG",])>0])]
intersect_names=intersect(names(col_name_BSG),names(col_name_TMPRSS2))
intersect_cells=table(Seurat::Idents(run.combined)[intersect_names])
count_matrix[15,names(intersect_cells)]=intersect_cells
TMPRSS2_cells=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["TMPRSS2",])>0])])
BSG_cells=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["BSG",])>0])])
count_matrix[16,names(BSG_cells)]=(intersect_cells/BSG_cells)*100
count_matrix[17,names(TMPRSS2_cells)]=(intersect_cells/TMPRSS2_cells)*100
new_data=data.frame("cells_percenrage"=c(((intersect_cells/TMPRSS2_cells)*100),((intersect_cells/BSG_cells)*100)),
genes=c(rep("TMPRSS2",length(TMPRSS2_cells)),rep("BSG",length(BSG_cells))),
cell_types=c(names(TMPRSS2_cells),names(BSG_cells)))
ggplot(data=new_data, aes(x=cell_types, y=cells_percenrage, fill=genes)) +
geom_bar(stat="identity") + theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

### Bar plot of intersect of both BSG and CTSL
col_name_CTSL=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["CTSL",])>0])]
col_name_BSG=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["BSG",])>0])]
intersect_names=intersect(names(col_name_BSG),names(col_name_CTSL))
intersect_cells=table(Seurat::Idents(run.combined)[intersect_names])
count_matrix[18,names(intersect_cells)]=intersect_cells
CTSL_cells=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["CTSL",])>0])])
BSG_cells=table(Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["BSG",])>0])])
count_matrix[19,names(BSG_cells)]=(intersect_cells/BSG_cells)*100
count_matrix[20,names(CTSL_cells)]=(intersect_cells/CTSL_cells)*100
new_data=data.frame("cells_percenrage"=c(((intersect_cells/CTSL_cells)*100),((intersect_cells/BSG_cells)*100)),
genes=c(rep("CTSL",length(CTSL_cells)),rep("BSG",length(BSG_cells))),
cell_types=c(names(CTSL_cells),names(BSG_cells)))
ggplot(data=new_data, aes(x=cell_types, y=cells_percenrage, fill=genes)) +
geom_bar(stat="identity") + theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

### Intersect of ACE2, TMPRSS2 and CTSL
col_name_TMPRSS2=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["TMPRSS2",])>0])]
col_name_ACE2=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["ACE2",])>0])]
col_name_CTSL=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["CTSL",])>0])]
intersect_names=intersect(names(col_name_ACE2),names(col_name_TMPRSS2))
intersect_names=intersect(intersect_names,names(col_name_CTSL))
intersect_cells=table(Seurat::Idents(run.combined)[intersect_names])
count_matrix[21,names(intersect_cells)]=intersect_cells
### Intersect of ACE2, TMPRSS2 and CTSB
col_name_TMPRSS2=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["TMPRSS2",])>0])]
col_name_ACE2=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["ACE2",])>0])]
col_name_CTSB=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["CTSB",])>0])]
intersect_names=intersect(names(col_name_ACE2),names(col_name_TMPRSS2))
intersect_names=intersect(intersect_names,names(col_name_CTSB))
intersect_cells=table(Seurat::Idents(run.combined)[intersect_names])
count_matrix[22,names(intersect_cells)]=intersect_cells
### Intersect of BSG, TMPRSS2 and CTSL
col_name_TMPRSS2=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["TMPRSS2",])>0])]
col_name_BSG=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["BSG",])>0])]
col_name_CTSL=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["CTSL",])>0])]
intersect_names=intersect(names(col_name_BSG),names(col_name_TMPRSS2))
intersect_names=intersect(intersect_names,names(col_name_CTSL))
intersect_cells=table(Seurat::Idents(run.combined)[intersect_names])
count_matrix[23,names(intersect_cells)]=intersect_cells
### Intersect of BSG, TMPRSS2 and CTSB
col_name_TMPRSS2=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["TMPRSS2",])>0])]
col_name_BSG=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["BSG",])>0])]
col_name_CTSB=Seurat::Idents(run.combined)[colnames(run.combined[,as.matrix(run.combined@assays$RNA["CTSB",])>0])]
intersect_names=intersect(names(col_name_BSG),names(col_name_TMPRSS2))
intersect_names=intersect(intersect_names,names(col_name_CTSB))
intersect_cells=table(Seurat::Idents(run.combined)[intersect_names])
count_matrix[24,names(intersect_cells)]=intersect_cells
count_matrix
## Immature neurons Sustentacular cells
## ACE2_count 0.0000000 15.000000
## ACE2_percent 0.0000000 46.875000
## TMPRSS2_count 5.0000000 345.000000
## TMPRSS2_percent 0.7507508 51.801802
## CTSL_count 167.0000000 360.000000
## CTSL_percent 10.2391171 22.072348
## BSG_count 347.0000000 622.000000
## BSG_percent 11.5129396 20.637027
## intersect_of_ACE2_TMPRSS2_count 0.0000000 9.000000
## ACE2_contri_percent NaN 60.000000
## TMPRSS2_contri_percent 0.0000000 2.608696
## intersect_of_ACE2_CTSL_count 0.0000000 11.000000
## ACE2_contri_percent NaN 73.333333
## CTSL_contri_percent 0.0000000 3.055556
## intersect_of_BSG_TMPRSS2_count 5.0000000 334.000000
## BSG_contri_percent 1.4409222 53.697749
## TMPRSS2_contri_percent 100.0000000 96.811594
## intersect_of_BSG_CTSL_count 148.0000000 347.000000
## BSG_contri_percent 42.6512968 55.787781
## CSTL_contri_percent 88.6227545 96.388889
## ACE2_TMPRSS2_CTSL 0.0000000 7.000000
## ACE2_TMPRSS2_CTSB 0.0000000 8.000000
## BSB_TMPRSS2_CTSL 2.0000000 217.000000
## BSG_TMPRSS2_CTSB 3.0000000 310.000000
## Olfactory HBCs Mature neurons
## ACE2_count 11.0000000 0.0000000
## ACE2_percent 34.3750000 0.0000000
## TMPRSS2_count 14.0000000 2.0000000
## TMPRSS2_percent 2.1021021 0.3003003
## CTSL_count 507.0000000 112.0000000
## CTSL_percent 31.0852238 6.8669528
## BSG_count 687.0000000 249.0000000
## BSG_percent 22.7936297 8.2614466
## intersect_of_ACE2_TMPRSS2_count 0.0000000 0.0000000
## ACE2_contri_percent 0.0000000 NaN
## TMPRSS2_contri_percent 0.0000000 0.0000000
## intersect_of_ACE2_CTSL_count 4.0000000 0.0000000
## ACE2_contri_percent 36.3636364 NaN
## CTSL_contri_percent 0.7889546 0.0000000
## intersect_of_BSG_TMPRSS2_count 13.0000000 2.0000000
## BSG_contri_percent 1.8922853 0.8032129
## TMPRSS2_contri_percent 92.8571429 100.0000000
## intersect_of_BSG_CTSL_count 342.0000000 109.0000000
## BSG_contri_percent 49.7816594 43.7751004
## CSTL_contri_percent 67.4556213 97.3214286
## ACE2_TMPRSS2_CTSL 0.0000000 0.0000000
## ACE2_TMPRSS2_CTSB 0.0000000 0.0000000
## BSB_TMPRSS2_CTSL 7.0000000 2.0000000
## BSG_TMPRSS2_CTSB 6.0000000 1.0000000
## Olfactory ensheathing glia Bowman's gland
## ACE2_count 0.0000000 4.0000000
## ACE2_percent 0.0000000 12.5000000
## TMPRSS2_count 6.0000000 239.0000000
## TMPRSS2_percent 0.9009009 35.8858859
## CTSL_count 90.0000000 281.0000000
## CTSL_percent 5.5180871 17.2286941
## BSG_count 157.0000000 772.0000000
## BSG_percent 5.2090246 25.6138023
## intersect_of_ACE2_TMPRSS2_count 0.0000000 1.0000000
## ACE2_contri_percent NaN 25.0000000
## TMPRSS2_contri_percent 0.0000000 0.4184100
## intersect_of_ACE2_CTSL_count 0.0000000 2.0000000
## ACE2_contri_percent NaN 50.0000000
## CTSL_contri_percent 0.0000000 0.7117438
## intersect_of_BSG_TMPRSS2_count 5.0000000 228.0000000
## BSG_contri_percent 3.1847134 29.5336788
## TMPRSS2_contri_percent 83.3333333 95.3974895
## intersect_of_BSG_CTSL_count 68.0000000 266.0000000
## BSG_contri_percent 43.3121019 34.4559585
## CSTL_contri_percent 75.5555556 94.6619217
## ACE2_TMPRSS2_CTSL 0.0000000 0.0000000
## ACE2_TMPRSS2_CTSB 0.0000000 1.0000000
## BSB_TMPRSS2_CTSL 2.0000000 91.0000000
## BSG_TMPRSS2_CTSB 1.0000000 215.0000000
## Olfactory microvillar cells GBCs
## ACE2_count 0.000000 2.000000
## ACE2_percent 0.000000 6.250000
## TMPRSS2_count 16.000000 39.000000
## TMPRSS2_percent 2.402402 5.855856
## CTSL_count 59.000000 55.000000
## CTSL_percent 3.617413 3.372164
## BSG_count 113.000000 67.000000
## BSG_percent 3.749171 2.222960
## intersect_of_ACE2_TMPRSS2_count 0.000000 2.000000
## ACE2_contri_percent NaN 100.000000
## TMPRSS2_contri_percent 0.000000 5.128205
## intersect_of_ACE2_CTSL_count 0.000000 1.000000
## ACE2_contri_percent NaN 50.000000
## CTSL_contri_percent 0.000000 1.818182
## intersect_of_BSG_TMPRSS2_count 16.000000 37.000000
## BSG_contri_percent 14.159292 55.223881
## TMPRSS2_contri_percent 100.000000 94.871795
## intersect_of_BSG_CTSL_count 57.000000 50.000000
## BSG_contri_percent 50.442478 74.626866
## CSTL_contri_percent 96.610169 90.909091
## ACE2_TMPRSS2_CTSL 0.000000 1.000000
## ACE2_TMPRSS2_CTSB 0.000000 1.000000
## BSB_TMPRSS2_CTSL 9.000000 32.000000
## BSG_TMPRSS2_CTSB 11.000000 33.000000
## Bar plots to show + expressed cells of some gene markers as shown beelow
### Bar side plot of ACE2, TMPRSS2 and (ACE2 and TMPRSS2)
count=c(count_matrix["ACE2_count",],count_matrix["TMPRSS2_count",],count_matrix["intersect_of_ACE2_TMPRSS2_count",])
count[count==0]<-1
A_T_AT_df=data.frame("count"=count,"cell_markers"= names(count),
"gene_markers"=c(rep("ACE2",8),rep("TMPRSS2",8),rep("Intersect (ACE2 and TMPRSS2)",8)))
ggplot(data=A_T_AT_df, aes(x=cell_markers, y=count, fill=gene_markers,width=.7)) +
geom_bar(stat="identity", position=position_dodge())+
theme_minimal() + scale_fill_manual(values=c('red','green',"blue")) +
scale_y_continuous(trans='log10') +
theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

### Bar side plot of BSG, CTSL and (BSG and CTSL)
count=c(count_matrix["BSG_count",],count_matrix["CTSL_count",],count_matrix["intersect_of_BSG_CTSL_count",])
count[count==0]<-1
B_C_BC_df=data.frame("count"=count,"cell_markers"= names(count),
"gene_markers"=c(rep("BSG",8),rep("CTSL",8),rep("Intersect (BSG and CTSL)",8)))
ggplot(data=B_C_BC_df, aes(x=cell_markers, y=count, fill=gene_markers,width=.7)) +
geom_bar(stat="identity", position=position_dodge())+
theme_minimal() + scale_fill_manual(values=c('red','green',"blue")) +
scale_y_continuous(trans='log10') +
theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

### Bar side plot of BSG, CTSL and (BSG and CTSL)
count=c(count_matrix["ACE2_count",],count_matrix["CTSL_count",],count_matrix["intersect_of_ACE2_CTSL_count",])
count[count==0]<-1
A_C_AC_df=data.frame("count"=count,"cell_markers"= names(count),
"gene_markers"=c(rep("ACE2",8),rep("CTSL",8),rep("Intersect (ACE2 and CTSL)",8)))
ggplot(data=A_C_AC_df, aes(x=cell_markers, y=count, fill=gene_markers,width=.7)) +
geom_bar(stat="identity", position=position_dodge())+
theme_minimal() + scale_fill_manual(values=c('red','green',"blue")) +
scale_y_continuous(trans='log10') +
theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

### Bar side plot of BSG, TMPRSS2 and (BSG and TMPRSS2)
count=c(count_matrix["BSG_count",],count_matrix["TMPRSS2_count",],count_matrix["intersect_of_BSG_TMPRSS2_count",])
count[count==0]<-1
B_T_BT_df=data.frame("count"=count,"cell_markers"= names(count),
"gene_markers"=c(rep("BSG",8),rep("TMPRSS2",8),rep("Intersect (BSG and TMPRSS2)",8)))
ggplot(data=B_T_BT_df, aes(x=cell_markers, y=count, fill=gene_markers,width=.7)) +
geom_bar(stat="identity", position=position_dodge())+
theme_minimal() + scale_fill_manual(values=c('red','green',"blue")) +
scale_y_continuous(trans='log10') +
theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

### Bar stack plot of ACE2, TMPRSS2 and (ACE2 and TMPRSS2)
count=c(count_matrix["ACE2_count",],count_matrix["TMPRSS2_count",],count_matrix["intersect_of_ACE2_TMPRSS2_count",])
count[count==0]<-1
A_T_AT_df=data.frame("count"=count,"cell_markers"= names(count),
"gene_markers"=c(rep("ACE2",8),rep("TMPRSS2",8),rep("Intersect (ACE2 and TMPRSS2)",8)))
ggplot(data=A_T_AT_df, aes(x=cell_markers, y=count, fill=gene_markers,width=.7)) +
geom_bar(stat="identity")+
theme_minimal() + scale_fill_manual(values=c('red','green',"blue")) +
scale_y_continuous(trans='log10') +
theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

### Bar stack plot of BSG, CTSL and (BSG and CTSL)
count=c(count_matrix["BSG_count",],count_matrix["CTSL_count",],count_matrix["intersect_of_BSG_CTSL_count",])
count[count==0]<-1
B_C_BC_df=data.frame("count"=count,"cell_markers"= names(count),
"gene_markers"=c(rep("BSG",8),rep("CTSL",8),rep("Intersect (BSG and CTSL)",8)))
ggplot(data=B_C_BC_df, aes(x=cell_markers, y=count, fill=gene_markers,width=.7)) +
geom_bar(stat="identity")+
theme_minimal() + scale_fill_manual(values=c('red','green',"blue")) +
scale_y_continuous(trans='log10') +
theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

### Bar stack plot of BSG, CTSL and (BSG and CTSL)
count=c(count_matrix["ACE2_count",],count_matrix["CTSL_count",],count_matrix["intersect_of_ACE2_CTSL_count",])
count[count==0]<-1
A_C_AC_df=data.frame("count"=count,"cell_markers"= names(count),
"gene_markers"=c(rep("ACE2",8),rep("CTSL",8),rep("Intersect (ACE2 and CTSL)",8)))
ggplot(data=A_C_AC_df, aes(x=cell_markers, y=count, fill=gene_markers,width=.7)) +
geom_bar(stat="identity")+
theme_minimal() + scale_fill_manual(values=c('red','green',"blue")) +
scale_y_continuous(trans='log10') +
theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

### Bar stack plot of BSG, TMPRSS2 and (BSG and TMPRSS2)
count=c(count_matrix["BSG_count",],count_matrix["TMPRSS2_count",],count_matrix["intersect_of_BSG_TMPRSS2_count",])
count[count==0]<-1
B_T_BT_df=data.frame("count"=count,"cell_markers"= names(count),
"gene_markers"=c(rep("BSG",8),rep("TMPRSS2",8),rep("Intersect (BSG and TMPRSS2)",8)))
ggplot(data=B_T_BT_df, aes(x=cell_markers, y=count, fill=gene_markers,width=.7)) +
geom_bar(stat="identity")+
theme_minimal() + scale_fill_manual(values=c('red','green',"blue")) +
scale_y_continuous(trans='log10') +
theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))

## DE genes analysis for selected cell types
run.combined=sce_3_1_1_after_5000_23
DefaultAssay(run.combined) <- "RNA"
seurat_subset_list=list()
p_val=0.05
l_fold=1
### DE genes analysis on ACE2+ and ACE2-
gene_marker="ACE2"
print(sort(count_matrix[1,]))
## Immature neurons Mature neurons
## 0 0
## Olfactory ensheathing glia Olfactory microvillar cells
## 0 0
## GBCs Bowman's gland
## 2 4
## Olfactory HBCs Sustentacular cells
## 11 15
selected_cells_type=sort(count_matrix[1,])
selected_cells_type=selected_cells_type[selected_cells_type>2]
de_methods=c("wilcox","MAST","bimod","t","negbinom","poisson","LR")
for(k in de_methods)
for(i in 1:length(selected_cells_type))
{
print(paste(k,"_",i,sep=""))
seurat_subset=subset(run.combined, idents = names(selected_cells_type)[i])
positive_cells_name = colnames(seurat_subset[,(as.matrix(seurat_subset@assays$RNA[gene_marker,])>0)])
negative_cells_name = colnames(seurat_subset[,(as.matrix(seurat_subset@assays$RNA[gene_marker,])<=0)])
pos_cells=paste(gene_marker,"+",sep="")
neg_cells=paste(gene_marker,"-",sep="")
pos=c(rep(pos_cells,length(positive_cells_name)))
names(pos)=positive_cells_name
neg=c(rep(neg_cells,length(negative_cells_name)))
names(neg)=negative_cells_name
pos_neg=c(pos,neg)
print(paste("cells in seurat object after subsetting with ",names(selected_cells_type)[i]," = ",dim(seurat_subset)[2],sep=""))
print(paste(names(selected_cells_type)[i]," total cells = ",sum(pos_neg[which(names(pos_neg)%in%rownames(seurat_subset@meta.data))]==pos_neg),sep=""))
print(paste(pos_cells," cells = ",sum(pos_neg==pos_cells),sep=""))
print(paste(neg_cells," cells = ",sum(pos_neg==neg_cells),sep=""))
seurat_subset=AddMetaData(
object = seurat_subset,
metadata = pos_neg,
col.name = 'seurat_subset'
)
seurat_subset_list[[paste(names(selected_cells_type)[i],"_",gene_marker,sep="")]]<-seurat_subset
Idents(seurat_subset)<-pos_neg
seurat_subset_genes=FindMarkers(seurat_subset, ident.1=pos_cells,ident.2=neg_cells,test.use = k,min.cells.group=1)
print(seurat_subset_genes[1:2,])
seurat_subset_genes_up=rownames(seurat_subset_genes)[seurat_subset_genes$avg_logFC>l_fold & seurat_subset_genes$p_val<p_val]
seurat_subset_genes_down=rownames(seurat_subset_genes)[seurat_subset_genes$avg_logFC<(-l_fold) & seurat_subset_genes$p_val<p_val]
seurat_subset_genes_names=c(seurat_subset_genes_up,seurat_subset_genes_down)
print(paste(names(selected_cells_type)[i],"_",pos_cells," genes = ",length(seurat_subset_genes_up),
" and ",names(selected_cells_type)[i],"_",neg_cells," genes = ",length(seurat_subset_genes_down),sep=""))
DE_info=data.frame("log2FC"=as.numeric(seurat_subset_genes$avg_logFC),"p_val"=as.numeric(seurat_subset_genes$p_val))
lab=rownames(seurat_subset_genes)
print(EnhancedVolcano(DE_info,
lab = lab,
x = 'log2FC',
y = 'p_val',
xlim = c(-2, 2),
title = names(selected_cells_type)[i],
pCutoff = p_val,
FCcutoff = l_fold,
colAlpha = 1))
}
## [1] "wilcox_1"
## [1] "cells in seurat object after subsetting with Bowman's gland = 903"
## [1] "Bowman's gland total cells = 903"
## [1] "ACE2+ cells = 4"
## [1] "ACE2- cells = 899"
## p_val avg_logFC pct.1 pct.2 p_val_adj
## AP003086.2 1.069865e-50 0.3138844 0.25 0.000 3.588112e-46
## AC009163.1 4.596423e-26 0.3134687 0.25 0.001 1.541548e-21
## [1] "Bowman's gland_ACE2+ genes = 3 and Bowman's gland_ACE2- genes = 2"

## [1] "wilcox_2"
## [1] "cells in seurat object after subsetting with Olfactory HBCs = 1087"
## [1] "Olfactory HBCs total cells = 1087"
## [1] "ACE2+ cells = 11"
## [1] "ACE2- cells = 1076"
## p_val avg_logFC pct.1 pct.2 p_val_adj
## SCN11A 1.905727e-15 0.3122231 0.182 0.004 6.391426e-11
## CHODL 1.521107e-07 0.2897168 0.182 0.010 5.101488e-03
## [1] "Olfactory HBCs_ACE2+ genes = 0 and Olfactory HBCs_ACE2- genes = 2"

## [1] "wilcox_3"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2+ cells = 15"
## [1] "ACE2- cells = 722"
## p_val avg_logFC pct.1 pct.2 p_val_adj
## CARTPT 1.257096e-05 0.2589717 0.2 0.022 0.4216048
## MYCL 2.385754e-05 0.6544713 0.8 0.303 0.8001343
## [1] "Sustentacular cells_ACE2+ genes = 0 and Sustentacular cells_ACE2- genes = 0"
## [1] "MAST_1"
## [1] "cells in seurat object after subsetting with Bowman's gland = 903"
## [1] "Bowman's gland total cells = 903"
## [1] "ACE2+ cells = 4"
## [1] "ACE2- cells = 899"
## Assuming data assay in position 1, with name et is log-transformed.
##
## Done!
## Combining coefficients and standard errors
## Warning in melt(coefAndCI, as.is = TRUE): The melt generic in data.table has
## been passed a array and will attempt to redirect to the relevant reshape2
## method; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(coefAndCI). In the next version, this warning will become an
## error.
## Calculating log-fold changes
## Warning in melt(lfc): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(lfc). In the
## next version, this warning will become an error.
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Warning in melt(llrt): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(llrt). In the
## next version, this warning will become an error.

## p_val avg_logFC pct.1 pct.2 p_val_adj
## R3HDM4 0.0007098688 0.7573622 1.0 0.205 1
## LINC01752 0.0008895931 0.4610704 0.5 0.007 1
## [1] "Bowman's gland_ACE2+ genes = 5 and Bowman's gland_ACE2- genes = 2"
## [1] "MAST_2"
## [1] "cells in seurat object after subsetting with Olfactory HBCs = 1087"
## [1] "Olfactory HBCs total cells = 1087"
## [1] "ACE2+ cells = 11"
## [1] "ACE2- cells = 1076"
## Assuming data assay in position 1, with name et is log-transformed.
##
## Done!
## Combining coefficients and standard errors
## Warning in melt(coefAndCI, as.is = TRUE): The melt generic in data.table has
## been passed a array and will attempt to redirect to the relevant reshape2
## method; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(coefAndCI). In the next version, this warning will become an
## error.
## Calculating log-fold changes
## Warning in melt(lfc): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(lfc). In the
## next version, this warning will become an error.
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Warning in melt(llrt): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(llrt). In the
## next version, this warning will become an error.

## p_val avg_logFC pct.1 pct.2 p_val_adj
## CTSA 9.283015e-06 0.4894748 0.182 0.243 0.3113337
## SMARCC2 3.723033e-05 0.5442944 0.273 0.235 1.0000000
## [1] "Olfactory HBCs_ACE2+ genes = 2 and Olfactory HBCs_ACE2- genes = 2"
## [1] "MAST_3"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2+ cells = 15"
## [1] "ACE2- cells = 722"
## Assuming data assay in position 1, with name et is log-transformed.
##
## Done!
## Combining coefficients and standard errors
## Warning in melt(coefAndCI, as.is = TRUE): The melt generic in data.table has
## been passed a array and will attempt to redirect to the relevant reshape2
## method; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(coefAndCI). In the next version, this warning will become an
## error.
## Calculating log-fold changes
## Warning in melt(lfc): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(lfc). In the
## next version, this warning will become an error.
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Warning in melt(llrt): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(llrt). In the
## next version, this warning will become an error.

## p_val avg_logFC pct.1 pct.2 p_val_adj
## AC009113.1 3.740250e-11 0.5956145 0.067 0.104 1.254405e-06
## C8orf58 4.486527e-05 0.5854481 0.200 0.112 1.000000e+00
## [1] "Sustentacular cells_ACE2+ genes = 1 and Sustentacular cells_ACE2- genes = 0"

## [1] "bimod_1"
## [1] "cells in seurat object after subsetting with Bowman's gland = 903"
## [1] "Bowman's gland total cells = 903"
## [1] "ACE2+ cells = 4"
## [1] "ACE2- cells = 899"
## p_val avg_logFC pct.1 pct.2 p_val_adj
## C19orf70 1.182295e-07 0.5201403 0.75 0.621 0.003965182
## PPP1R32 6.583710e-07 0.3598497 0.50 0.012 0.022080447
## [1] "Bowman's gland_ACE2+ genes = 2 and Bowman's gland_ACE2- genes = 3"

## [1] "bimod_2"
## [1] "cells in seurat object after subsetting with Olfactory HBCs = 1087"
## [1] "Olfactory HBCs total cells = 1087"
## [1] "ACE2+ cells = 11"
## [1] "ACE2- cells = 1076"
## p_val avg_logFC pct.1 pct.2 p_val_adj
## IGLC3 6.825386e-13 2.190882 0.364 0.205 2.289098e-08
## IGLC2 5.545072e-10 1.824826 0.364 0.398 1.859706e-05
## [1] "Olfactory HBCs_ACE2+ genes = 2 and Olfactory HBCs_ACE2- genes = 2"

## [1] "bimod_3"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2+ cells = 15"
## [1] "ACE2- cells = 722"
## p_val avg_logFC pct.1 pct.2 p_val_adj
## URB1-AS1 2.788252e-16 0.8820008 0.600 0.269 9.351240e-12
## AL034549.1 2.358246e-15 0.6249744 0.267 0.118 7.909085e-11
## [1] "Sustentacular cells_ACE2+ genes = 1 and Sustentacular cells_ACE2- genes = 2"

## [1] "t_1"
## [1] "cells in seurat object after subsetting with Bowman's gland = 903"
## [1] "Bowman's gland total cells = 903"
## [1] "ACE2+ cells = 4"
## [1] "ACE2- cells = 899"
## p_val avg_logFC pct.1 pct.2 p_val_adj
## HDLBP 1.724021e-155 -0.8721947 0 0.638 5.782021e-151
## MPHOSPH8 1.049661e-153 -0.8076480 0 0.628 3.520352e-149
## [1] "Bowman's gland_ACE2+ genes = 0 and Bowman's gland_ACE2- genes = 9"

## [1] "t_2"
## [1] "cells in seurat object after subsetting with Olfactory HBCs = 1087"
## [1] "Olfactory HBCs total cells = 1087"
## [1] "ACE2+ cells = 11"
## [1] "ACE2- cells = 1076"
## p_val avg_logFC pct.1 pct.2 p_val_adj
## LAMB1 9.368870e-135 -0.7778893 0 0.486 3.142132e-130
## ATP13A3 1.615811e-103 -0.5378039 0 0.393 5.419108e-99
## [1] "Olfactory HBCs_ACE2+ genes = 0 and Olfactory HBCs_ACE2- genes = 3"

## [1] "t_3"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2+ cells = 15"
## [1] "ACE2- cells = 722"
## p_val avg_logFC pct.1 pct.2 p_val_adj
## TMEM238 8.163443e-25 -0.3774210 0 0.190 2.737855e-20
## XIST 1.701171e-17 -0.5124281 0 0.109 5.705389e-13
## [1] "Sustentacular cells_ACE2+ genes = 0 and Sustentacular cells_ACE2- genes = 2"
## [1] "negbinom_1"
## [1] "cells in seurat object after subsetting with Bowman's gland = 903"
## [1] "Bowman's gland total cells = 903"
## [1] "ACE2+ cells = 4"
## [1] "ACE2- cells = 899"
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in FUN(X[[i]], ...): Skipping gene --- 649. Fewer than 3 cells in both
## clusters.
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## p_val avg_logFC pct.1 pct.2 p_val_adj
## LINC01752 1.263145e-07 0.3988132 0.5 0.007 0.004236336
## PPP1R32 1.402014e-06 0.3933035 0.5 0.012 0.047020752
## [1] "Bowman's gland_ACE2+ genes = 0 and Bowman's gland_ACE2- genes = 7"
## [1] "negbinom_2"
## [1] "cells in seurat object after subsetting with Olfactory HBCs = 1087"
## [1] "Olfactory HBCs total cells = 1087"
## [1] "ACE2+ cells = 11"
## [1] "ACE2- cells = 1076"
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in glm.nb(formula = fmla, data = latent.vars): alternation limit reached

## p_val avg_logFC pct.1 pct.2 p_val_adj
## IGLC3 1.960891e-10 1.797373 0.364 0.205 6.576436e-06
## MT-ND1 9.184367e-05 -1.105728 1.000 0.995 1.000000e+00
## [1] "Olfactory HBCs_ACE2+ genes = 2 and Olfactory HBCs_ACE2- genes = 9"
## [1] "negbinom_3"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2+ cells = 15"
## [1] "ACE2- cells = 722"
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in glm.nb(formula = fmla, data = latent.vars): alternation limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## p_val avg_logFC pct.1 pct.2 p_val_adj
## BICDL1 6.437705e-08 0.3978762 0.467 0.071 0.002159077
## MRPL48 9.001799e-07 0.6287818 0.733 0.355 0.030190234
## [1] "Sustentacular cells_ACE2+ genes = 5 and Sustentacular cells_ACE2- genes = 2"
## [1] "poisson_1"
## [1] "cells in seurat object after subsetting with Bowman's gland = 903"
## [1] "Bowman's gland total cells = 903"
## [1] "ACE2+ cells = 4"
## [1] "ACE2- cells = 899"
## Warning in FUN(X[[i]], ...): Skipping gene --- 649. Fewer than 3 cells in both
## clusters.

## p_val avg_logFC pct.1 pct.2 p_val_adj
## SAA1 7.245093e-219 2.152606 0.25 0.271 2.429859e-214
## LYZ 3.583759e-134 -3.023938 0.75 0.934 1.201921e-129
## [1] "Bowman's gland_ACE2+ genes = 14 and Bowman's gland_ACE2- genes = 23"

## [1] "poisson_2"
## [1] "cells in seurat object after subsetting with Olfactory HBCs = 1087"
## [1] "Olfactory HBCs total cells = 1087"
## [1] "ACE2+ cells = 11"
## [1] "ACE2- cells = 1076"
## p_val avg_logFC pct.1 pct.2 p_val_adj
## IGLC3 1.832406e-132 1.797373 0.364 0.205 6.145523e-128
## IGLC2 1.357811e-65 1.425866 0.364 0.398 4.553827e-61
## [1] "Olfactory HBCs_ACE2+ genes = 2 and Olfactory HBCs_ACE2- genes = 9"
## [1] "poisson_3"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2+ cells = 15"
## [1] "ACE2- cells = 722"
## p_val avg_logFC pct.1 pct.2 p_val_adj
## BPIFB1 0.000000e+00 2.5050723 0.133 0.062 0.000000e+00
## GSTP1 4.435057e-238 0.4987906 1.000 0.992 1.487429e-233
## [1] "Sustentacular cells_ACE2+ genes = 14 and Sustentacular cells_ACE2- genes = 2"
## Warning: One or more p-values is 0. Converting to 10^-1 * current lowest non-
## zero p-value...

## [1] "LR_1"
## [1] "cells in seurat object after subsetting with Bowman's gland = 903"
## [1] "Bowman's gland total cells = 903"
## [1] "ACE2+ cells = 4"
## [1] "ACE2- cells = 899"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## p_val avg_logFC pct.1 pct.2 p_val_adj
## LINC01752 6.320178e-05 0.4610704 0.50 0.007 1
## C3orf33 3.938426e-04 0.5344172 0.75 0.053 1
## [1] "Bowman's gland_ACE2+ genes = 4 and Bowman's gland_ACE2- genes = 3"
## [1] "LR_2"
## [1] "cells in seurat object after subsetting with Olfactory HBCs = 1087"
## [1] "Olfactory HBCs total cells = 1087"
## [1] "ACE2+ cells = 11"
## [1] "ACE2- cells = 1076"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## p_val avg_logFC pct.1 pct.2 p_val_adj
## TULP3 8.050212e-05 0.7332917 0.455 0.096 1
## ARL4D 1.345051e-04 -1.1427918 0.091 0.589 1
## [1] "Olfactory HBCs_ACE2+ genes = 1 and Olfactory HBCs_ACE2- genes = 3"
## [1] "LR_3"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2+ cells = 15"
## [1] "ACE2- cells = 722"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## p_val avg_logFC pct.1 pct.2 p_val_adj
## MYCL 0.0004320384 0.6544713 0.8 0.303 1
## URB1-AS1 0.0004526378 0.8820008 0.6 0.269 1
## [1] "Sustentacular cells_ACE2+ genes = 1 and Sustentacular cells_ACE2- genes = 1"

### DE genes analysis on intersection of (ACE2 and TMPRSS2)+ and (ACE2 and TMPRSS2)-
gene_marker=c("ACE2","TMPRSS2")
print(sort(count_matrix[1,]))
## Immature neurons Mature neurons
## 0 0
## Olfactory ensheathing glia Olfactory microvillar cells
## 0 0
## GBCs Bowman's gland
## 2 4
## Olfactory HBCs Sustentacular cells
## 11 15
selected_cells_type=sort(count_matrix[9,])
selected_cells_type=selected_cells_type[selected_cells_type>2]
for(k in de_methods)
for(i in 1:length(selected_cells_type))
{
print(paste(k,"_",i,sep=""))
seurat_subset=subset(run.combined, idents = names(selected_cells_type)[i])
positive_cells_name = colnames(seurat_subset[,(as.matrix(seurat_subset@assays$RNA[gene_marker[1],])>0 & as.matrix(seurat_subset@assays$RNA[gene_marker[2],])>0)])
negative_cells_name = colnames(seurat_subset[,(as.matrix(seurat_subset@assays$RNA[gene_marker[1],])<=0 | as.matrix(seurat_subset@assays$RNA[gene_marker[2],])<=0)])
pos_cells=paste(gene_marker[1],"_and_",gene_marker[2],"+",sep="")
neg_cells=paste(gene_marker[1],"_and_",gene_marker[2],"-",sep="")
pos=c(rep(pos_cells,length(positive_cells_name)))
names(pos)=positive_cells_name
neg=c(rep(neg_cells,length(negative_cells_name)))
names(neg)=negative_cells_name
pos_neg=c(pos,neg)
print(paste("cells in seurat object after subsetting with ",names(selected_cells_type)[i]," = ",dim(seurat_subset)[2],sep=""))
print(paste(names(selected_cells_type)[i]," total cells = ",sum(pos_neg[which(names(pos_neg)%in%rownames(seurat_subset@meta.data))]==pos_neg),sep=""))
print(paste(pos_cells," cells = ",sum(pos_neg==pos_cells),sep=""))
print(paste(neg_cells," cells = ",sum(pos_neg==neg_cells),sep=""))
seurat_subset=AddMetaData(
object = seurat_subset,
metadata = pos_neg,
col.name = 'seurat_subset'
)
seurat_subset_list[[paste(names(selected_cells_type)[i],"_",gene_marker[1],"_and_",gene_marker[2],sep="")]]<-seurat_subset
Idents(seurat_subset)<-pos_neg
seurat_subset_genes=FindMarkers(seurat_subset, ident.1=pos_cells,ident.2=neg_cells,test.use = k,min.cells.group=1)
print(seurat_subset_genes[1:2,])
seurat_subset_genes_up=rownames(seurat_subset_genes)[seurat_subset_genes$avg_logFC>l_fold & seurat_subset_genes$p_val<p_val]
seurat_subset_genes_down=rownames(seurat_subset_genes)[seurat_subset_genes$avg_logFC<(-l_fold) & seurat_subset_genes$p_val<p_val]
seurat_subset_genes_names=c(seurat_subset_genes_up,seurat_subset_genes_down)
print(paste(names(selected_cells_type)[i],"_",pos_cells," genes = ",length(seurat_subset_genes_up),
" and ",names(selected_cells_type)[i],"_",neg_cells," genes = ",length(seurat_subset_genes_down),sep=""))
DE_info=data.frame("log2FC"=as.numeric(seurat_subset_genes$avg_logFC),"p_val"=as.numeric(seurat_subset_genes$p_val))
lab=rownames(seurat_subset_genes)
print(EnhancedVolcano(DE_info,
lab = lab,
x = 'log2FC',
y = 'p_val',
xlim = c(-5, 5),
title = names(selected_cells_type)[i],
pCutoff = p_val,
FCcutoff = l_fold,
colAlpha = 1))
}
## [1] "wilcox_1"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2_and_TMPRSS2+ cells = 9"
## [1] "ACE2_and_TMPRSS2- cells = 728"
## p_val avg_logFC pct.1 pct.2 p_val_adj
## PLXNA1 1.104891e-06 0.2963436 0.778 0.169 0.03705582
## MTX3 2.878126e-06 0.3344437 0.778 0.201 0.09652658
## [1] "Sustentacular cells_ACE2_and_TMPRSS2+ genes = 0 and Sustentacular cells_ACE2_and_TMPRSS2- genes = 0"
## [1] "MAST_1"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2_and_TMPRSS2+ cells = 9"
## [1] "ACE2_and_TMPRSS2- cells = 728"
## Assuming data assay in position 1, with name et is log-transformed.
##
## Done!
## Combining coefficients and standard errors
## Warning in melt(coefAndCI, as.is = TRUE): The melt generic in data.table has
## been passed a array and will attempt to redirect to the relevant reshape2
## method; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(coefAndCI). In the next version, this warning will become an
## error.
## Calculating log-fold changes
## Warning in melt(lfc): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(lfc). In the
## next version, this warning will become an error.
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Warning in melt(llrt): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(llrt). In the
## next version, this warning will become an error.

## p_val avg_logFC pct.1 pct.2 p_val_adj
## C1orf159 1.818218e-09 0.7760622 0.111 0.081 6.097939e-05
## C8orf58 6.750627e-09 0.8492568 0.222 0.113 2.264025e-04
## [1] "Sustentacular cells_ACE2_and_TMPRSS2+ genes = 3 and Sustentacular cells_ACE2_and_TMPRSS2- genes = 0"

## [1] "bimod_1"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2_and_TMPRSS2+ cells = 9"
## [1] "ACE2_and_TMPRSS2- cells = 728"
## p_val avg_logFC pct.1 pct.2 p_val_adj
## SELENOO 5.270310e-15 0.7256710 0.333 0.209 1.767557e-10
## FAM19A4 2.976249e-12 0.7768982 0.222 0.098 9.981744e-08
## [1] "Sustentacular cells_ACE2_and_TMPRSS2+ genes = 5 and Sustentacular cells_ACE2_and_TMPRSS2- genes = 2"

## [1] "t_1"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2_and_TMPRSS2+ cells = 9"
## [1] "ACE2_and_TMPRSS2- cells = 728"
## p_val avg_logFC pct.1 pct.2 p_val_adj
## MYL9 9.828220e-26 -0.2544036 0 0.227 3.296189e-21
## TMEM238 8.437425e-25 -0.3748267 0 0.188 2.829743e-20
## [1] "Sustentacular cells_ACE2_and_TMPRSS2+ genes = 0 and Sustentacular cells_ACE2_and_TMPRSS2- genes = 2"
## [1] "negbinom_1"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2_and_TMPRSS2+ cells = 9"
## [1] "ACE2_and_TMPRSS2- cells = 728"
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in glm.nb(formula = fmla, data = latent.vars): alternation limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in FUN(X[[i]], ...): Skipping gene --- 1650. Fewer than 3 cells in both
## clusters.
## Warning in FUN(X[[i]], ...): Skipping gene --- 1780. Fewer than 3 cells in both
## clusters.
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in FUN(X[[i]], ...): Skipping gene --- 1957. Fewer than 3 cells in both
## clusters.
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## p_val avg_logFC pct.1 pct.2 p_val_adj
## KLHL32 4.679348e-07 0.2767530 0.333 0.011 0.01569360
## NBPF15 1.775853e-06 0.4278287 0.444 0.084 0.05955854
## [1] "Sustentacular cells_ACE2_and_TMPRSS2+ genes = 9 and Sustentacular cells_ACE2_and_TMPRSS2- genes = 0"
## [1] "poisson_1"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2_and_TMPRSS2+ cells = 9"
## [1] "ACE2_and_TMPRSS2- cells = 728"
## Warning in FUN(X[[i]], ...): Skipping gene --- 1650. Fewer than 3 cells in both
## clusters.
## Warning in FUN(X[[i]], ...): Skipping gene --- 1780. Fewer than 3 cells in both
## clusters.
## Warning in FUN(X[[i]], ...): Skipping gene --- 1957. Fewer than 3 cells in both
## clusters.
## p_val avg_logFC pct.1 pct.2 p_val_adj
## SCGB1A1 0 3.269476 0.111 0.058 0
## BPIFB1 0 2.985891 0.111 0.063 0
## [1] "Sustentacular cells_ACE2_and_TMPRSS2+ genes = 27 and Sustentacular cells_ACE2_and_TMPRSS2- genes = 2"
## Warning: One or more p-values is 0. Converting to 10^-1 * current lowest non-
## zero p-value...

## [1] "LR_1"
## [1] "cells in seurat object after subsetting with Sustentacular cells = 737"
## [1] "Sustentacular cells total cells = 737"
## [1] "ACE2_and_TMPRSS2+ cells = 9"
## [1] "ACE2_and_TMPRSS2- cells = 728"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## p_val avg_logFC pct.1 pct.2 p_val_adj
## COPS7A 0.0004655404 0.7991358 0.889 0.438 1
## CCDC184 0.0004961425 0.3207393 0.333 0.052 1
## [1] "Sustentacular cells_ACE2_and_TMPRSS2+ genes = 3 and Sustentacular cells_ACE2_and_TMPRSS2- genes = 0"

## Finding Stouffer score after log transfer and zscore
### Loading seaurat object
sce_3_1_1_after_5000_23 <- readRDS("~/corona_project/sce_3_1_1_after_5000_23.rds")
mat<-sce_3_1_1_after_5000_23@assays$RNA@counts
### Cell filtering
c_s <- Matrix::colSums(mat>=3)
l1 = stats::quantile(c_s,probs = 0.001)
l2 = stats::quantile(c_s,probs = 1)
keep_cells = intersect(which(c_s>=l1), which(c_s<=l2))
mat<-mat[,keep_cells]
### Gene filtering
r_s <- Matrix::rowSums(mat > 2)
keep_genes = which(r_s > 3)
mat<-mat[keep_genes,]
### Median normalization
r_s<-Matrix::rowSums(t(mat))
r_s_med<-stats::median(r_s)
mat<-Matrix::t(t(mat)/(r_s/r_s_med))
### Loading media genes file and intersecting genes and filtering matrix common genes
media_genes <- read_excel("~/corona_project/media-6.xlsx")$...3[-1]
## New names:
## * `` -> ...2
## * `` -> ...3
## * `` -> ...5
## * `` -> ...6
## * `` -> ...7
## * ...
common_genes=intersect(rownames(mat),media_genes)
mat=mat[common_genes,]
### Further insuring have cells with atleast 10% expressed genes
mat<-mat[,Matrix::colSums(mat>0)>(dim(mat)[1]/10)]
print(dim(mat))
## [1] 293 3493
seurat_RNA_mat<-mat
print(sum(apply(seurat_RNA_mat,1,function(x) sum(x)==0)))
## [1] 0
print(sum(apply(seurat_RNA_mat,2,function(x) sum(x)==0)))
## [1] 0
print(dim(seurat_RNA_mat))
## [1] 293 3493
### First psedo count 1 then log2 then zscore using scale function in R
seurat_RNA_mat<-scale(log2(seurat_RNA_mat+1),center = TRUE, scale = TRUE)
### Stouffer score using corto package
stouffer_score<- apply(seurat_RNA_mat,2,corto::stouffer)
print(sum(is.na(stouffer_score)))
## [1] 0
cell_types<-Idents(sce_3_1_1_after_5000_23)[names(stouffer_score)]
unique_cell_types=unique(cell_types)
Stouffer_score_df<-data.frame("Stouffer_score"=stouffer_score, "cell_types"=cell_types)
med=c()
med_n=c()
ct=unique(cell_types)
for(i in ct)
{
med = c(med,median(stouffer_score[cell_types==i]))
med_n = c(med_n,i)
}
med_df=data.frame("cell_types"=med_n, "med"=med)
ggplot(med_df, aes(x=cell_types, y=med, fill=cell_types)) +
geom_bar(stat="identity")+theme_minimal() +
theme_classic()+ theme(axis.text.x = element_text(angle = 90, hjust = 1))
